{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rank Features (Solution)\n",
    "\n",
    "The creator of Shapley Additive Explanations, Scott Lundberg, has written an efficient implementation that we can install and use.  We’ll be able to use this to determine both local feature importance (for a single observation) and global feature importance (for all training samples as a whole).  To aggregate local feature importance into global feature importance, we take the absolute values of the local feature importances, and then average them.\n",
    "\n",
    "We can calculate the feature importance using sklearn and using the Shap library.\n",
    "\n",
    "Based on the feature importances, we can think about modifying features to improve them.  Then we can re-train the model on the modified features.\n",
    "\n",
    "Finally, we can prune the feature set to just use the most relevant features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Note, this will install zipline and alphalens, which will take some time\n",
    "import sys\n",
    "!{sys.executable} -m pip install --quiet -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import time\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = (14, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import project_helper\n",
    "from zipline.data import bundles\n",
    "\n",
    "os.environ['ZIPLINE_ROOT'] = os.path.join(os.getcwd(), '..', '..', 'data', 'module_4_quizzes_eod')\n",
    "\n",
    "ingest_func = bundles.csvdir.csvdir_equities(['daily'], project_helper.EOD_BUNDLE_NAME)\n",
    "bundles.register(project_helper.EOD_BUNDLE_NAME, ingest_func)\n",
    "\n",
    "print('Data Registered')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from zipline.pipeline import Pipeline\n",
    "from zipline.pipeline.factors import AverageDollarVolume\n",
    "from zipline.utils.calendars import get_calendar\n",
    "\n",
    "\n",
    "universe = AverageDollarVolume(window_length=120).top(500) \n",
    "trading_calendar = get_calendar('NYSE') \n",
    "bundle_data = bundles.load(project_helper.EOD_BUNDLE_NAME)\n",
    "engine = project_helper.build_pipeline_engine(bundle_data, trading_calendar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test\n",
    "universe_end_date = pd.Timestamp('2016-01-05', tz='UTC')\n",
    "\n",
    "universe_tickers = engine\\\n",
    "    .run_pipeline(\n",
    "        Pipeline(screen=universe),\n",
    "        universe_end_date,\n",
    "        universe_end_date)\\\n",
    "    .index.get_level_values(1)\\\n",
    "    .values.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from zipline.data.data_portal import DataPortal\n",
    "\n",
    "data_portal = DataPortal(\n",
    "    bundle_data.asset_finder,\n",
    "    trading_calendar=trading_calendar,\n",
    "    first_trading_day=bundle_data.equity_daily_bar_reader.first_trading_day,\n",
    "    equity_minute_reader=None,\n",
    "    equity_daily_reader=bundle_data.equity_daily_bar_reader,\n",
    "    adjustment_reader=bundle_data.adjustment_reader)\n",
    "\n",
    "def get_pricing(data_portal, trading_calendar, assets, start_date, end_date, field='close'):\n",
    "    end_dt = pd.Timestamp(end_date.strftime('%Y-%m-%d'), tz='UTC', offset='C')\n",
    "    start_dt = pd.Timestamp(start_date.strftime('%Y-%m-%d'), tz='UTC', offset='C')\n",
    "\n",
    "    end_loc = trading_calendar.closes.index.get_loc(end_dt)\n",
    "    start_loc = trading_calendar.closes.index.get_loc(start_dt)\n",
    "\n",
    "    return data_portal.get_history_window(\n",
    "        assets=assets,\n",
    "        end_dt=end_dt,\n",
    "        bar_count=end_loc - start_loc,\n",
    "        frequency='1d',\n",
    "        field=field,\n",
    "        data_frequency='daily')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make Factors\n",
    "\n",
    "\n",
    "- Take the same factors we have been using:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from zipline.pipeline.factors import CustomFactor, DailyReturns, Returns, SimpleMovingAverage\n",
    "from zipline.pipeline.data import USEquityPricing\n",
    "\n",
    "factor_start_date = universe_end_date - pd.DateOffset(years=3, days=2)\n",
    "sector = project_helper.Sector()\n",
    "\n",
    "def momentum_1yr(window_length, universe, sector):\n",
    "    return Returns(window_length=window_length, mask=universe) \\\n",
    "        .demean(groupby=sector) \\\n",
    "        .rank() \\\n",
    "        .zscore()\n",
    "\n",
    "def mean_reversion_5day_sector_neutral(window_length, universe, sector):\n",
    "    return -Returns(window_length=window_length, mask=universe) \\\n",
    "        .demean(groupby=sector) \\\n",
    "        .rank() \\\n",
    "        .zscore()\n",
    "\n",
    "def mean_reversion_5day_sector_neutral_smoothed(window_length, universe, sector):\n",
    "    unsmoothed_factor = mean_reversion_5day_sector_neutral(window_length, universe, sector)\n",
    "    return SimpleMovingAverage(inputs=[unsmoothed_factor], window_length=window_length) \\\n",
    "        .rank() \\\n",
    "        .zscore()\n",
    "\n",
    "class CTO(Returns):\n",
    "    \"\"\"\n",
    "    Computes the overnight return, per hypothesis from\n",
    "    https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2554010\n",
    "    \"\"\"\n",
    "    inputs = [USEquityPricing.open, USEquityPricing.close]\n",
    "    \n",
    "    def compute(self, today, assets, out, opens, closes):\n",
    "        \"\"\"\n",
    "        The opens and closes matrix is 2 rows x N assets, with the most recent at the bottom.\n",
    "        As such, opens[-1] is the most recent open, and closes[0] is the earlier close\n",
    "        \"\"\"\n",
    "        out[:] = (opens[-1] - closes[0]) / closes[0]\n",
    "\n",
    "        \n",
    "class TrailingOvernightReturns(Returns):\n",
    "    \"\"\"\n",
    "    Sum of trailing 1m O/N returns\n",
    "    \"\"\"\n",
    "    window_safe = True\n",
    "    \n",
    "    def compute(self, today, asset_ids, out, cto):\n",
    "        out[:] = np.nansum(cto, axis=0)\n",
    "\n",
    "        \n",
    "def overnight_sentiment(cto_window_length, trail_overnight_returns_window_length, universe):\n",
    "    cto_out = CTO(mask=universe, window_length=cto_window_length)\n",
    "    return TrailingOvernightReturns(inputs=[cto_out], window_length=trail_overnight_returns_window_length) \\\n",
    "        .rank() \\\n",
    "        .zscore()\n",
    "\n",
    "def overnight_sentiment_smoothed(cto_window_length, trail_overnight_returns_window_length, universe):\n",
    "    unsmoothed_factor = overnight_sentiment(cto_window_length, trail_overnight_returns_window_length, universe)\n",
    "    return SimpleMovingAverage(inputs=[unsmoothed_factor], window_length=trail_overnight_returns_window_length) \\\n",
    "        .rank() \\\n",
    "        .zscore()\n",
    "\n",
    "universe = AverageDollarVolume(window_length=120).top(500)\n",
    "sector = project_helper.Sector()\n",
    "\n",
    "pipeline = Pipeline(screen=universe)\n",
    "pipeline.add(\n",
    "    momentum_1yr(252, universe, sector),\n",
    "    'Momentum_1YR')\n",
    "pipeline.add(\n",
    "    mean_reversion_5day_sector_neutral_smoothed(20, universe, sector),\n",
    "    'Mean_Reversion_Sector_Neutral_Smoothed')\n",
    "pipeline.add(\n",
    "    overnight_sentiment_smoothed(2, 10, universe),\n",
    "    'Overnight_Sentiment_Smoothed')\n",
    "\n",
    "all_factors = engine.run_pipeline(pipeline, factor_start_date, universe_end_date)\n",
    "\n",
    "all_factors.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add sector code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipeline.add(sector, 'sector_code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Universal Quant Features\n",
    "\n",
    "* stock volatility: zipline has a custom factor called AnnualizedVolatility.  The [source code is here](https://github.com/quantopian/zipline/blob/master/zipline/pipeline/factors/basic.py) and also pasted below:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Annualized volatility.\n",
    "Create `AnnualizedVolatility` objects for 20 day and 120 day (one month and six-month) time windows.  Remember to set the `mask` parameter to the `universe` object created earlier (this filters the stocks to match the list in the `universe`).  Convert these to ranks, and then convert the ranks to zscores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from zipline.pipeline.factors import AnnualizedVolatility\n",
    "volatility_20d = AnnualizedVolatility(window_length=20, mask=universe).rank().zscore()\n",
    "volatility_120d = AnnualizedVolatility(window_length=120, mask=universe).rank().zscore()\n",
    "pipeline.add(volatility_20d, 'volatility_20d')\n",
    "pipeline.add(volatility_120d, 'volatility_120d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Average Dollar Volume feature\n",
    "[AverageDollarVolume](http://www.zipline.io/appendix.html#zipline.pipeline.factors.AverageDollarVolume):\n",
    "Use 20 day and 120 day `window_length`, rank and then zscore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from zipline.pipeline.factors import AverageDollarVolume # already imported earlier, but shown here for reference\n",
    "adv_20d = AverageDollarVolume(window_length=20, mask=universe).rank().zscore()\n",
    "adv_120d = AverageDollarVolume(window_length=120, mask=universe).rank().zscore()\n",
    "pipeline.add(adv_20d, 'adv_20d')\n",
    "pipeline.add(adv_120d, 'adv_120d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regime Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### market dispersion feature\n",
    "\n",
    "Calculate the mean returns\n",
    "\n",
    "$\\mu = \\sum_{t=0}^{T}\\sum_{i=1}^{N}r_{i,t}$\n",
    "\n",
    "$\\sqrt{\\frac{1}{T} \\sum_{t=0}^{T}  \\frac{1}{N}\\sum_{i=1}^{N}(r_{i,t} - \\mu)^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MarketDispersion(CustomFactor):\n",
    "    inputs = [DailyReturns()]\n",
    "    window_length = 1\n",
    "    window_safe = True\n",
    "\n",
    "    def compute(self, today, assets, out, returns):\n",
    "        # returns are days in rows, assets across columns\n",
    "        mean_returns = np.nanmean(returns)\n",
    "        out[:] = np.sqrt(np.nanmean((returns - mean_returns)**2))\n",
    "        \n",
    "pipeline.add(SimpleMovingAverage(inputs=[MarketDispersion(mask=universe)], window_length=20), 'dispersion_20d')\n",
    "pipeline.add(SimpleMovingAverage(inputs=[MarketDispersion(mask=universe)], window_length=120), 'dispersion_120d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Market volatility feature\n",
    "* High and low volatility  \n",
    "We'll also build a class for market volatility, which inherits from [CustomFactor](http://www.zipline.io/appendix.html?highlight=customfactor#zipline.pipeline.CustomFactor).\n",
    "\n",
    "##### Market return\n",
    "$r_{m,t} = \\sum_{i=1}^{N}r_{i,t}$ for each day $t$ in `window_length`.  \n",
    "\n",
    "##### Average market return\n",
    "Also calculate the average market return over the `window_length` $W$ of days:  \n",
    "$\\mu_{m} = \\frac{1}{N}\\sum_{t=1}^{T} r_{m,t}$\n",
    "\n",
    "#### Standard deviation of market return\n",
    "Then calculate the standard deviation of the market return  \n",
    "$\\sigma_{m,t} = \\sqrt{252 \\times \\frac{1}{N} \\sum_{t=1}^{T}(r_{m,t} - \\mu_{m})^2 } $ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MarketVolatility(CustomFactor):\n",
    "    inputs = [DailyReturns()]\n",
    "    window_length = 1  # We'll want to set this in the constructor when creating the object.\n",
    "    window_safe = True\n",
    "    \n",
    "    \n",
    "    def compute(self, today, assets, out, returns):\n",
    "        DAILY_TO_ANNUAL_SCALAR = 252.  # 252 trading days in a year\n",
    "        \"\"\"\n",
    "        For each row (each row represents one day of returns), \n",
    "        calculate the average of the cross-section of stock returns\n",
    "        So that market_returns has one value for each day in the window_length\n",
    "        So choose the appropriate axis (please see hints above)\n",
    "        \"\"\"\n",
    "        mkt_returns = np.nanmean(returns, axis=1) \n",
    "        \n",
    "        \"\"\" \n",
    "        Calculate the mean of market returns\n",
    "        \"\"\"\n",
    "        mkt_returns_mu = np.nanmean(mkt_returns)\n",
    "        \n",
    "        \"\"\"\n",
    "        Calculate the standard deviation of the market returns, then annualize them.\n",
    "        \"\"\"\n",
    "        out[:] = np.sqrt(DAILY_TO_ANNUAL_SCALAR * np.nanmean((mkt_returns-mkt_returns_mu)**2))\n",
    "        \n",
    "# create market volatility features using one month and six-month windows\n",
    "market_vol_20d = MarketVolatility(window_length=20)\n",
    "market_vol_120d = MarketVolatility(window_length=120)\n",
    "\n",
    "# add market volatility features to pipeline\n",
    "pipeline.add(market_vol_20d, 'market_vol_20d')\n",
    "pipeline.add(market_vol_120d, 'market_vol_120d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run pipeline to calculate features\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors = engine.run_pipeline(pipeline, factor_start_date, universe_end_date)\n",
    "all_factors.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make Date Parts\n",
    "* we make colums to for the trees to split on that might capture trader/investor behavior due to calendar anomalies.\n",
    "* We can get the dates from the index of the dataframe that is returned from running the pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## January, December\n",
    "* Create a numpy array that has 1 when the month is January, and 0 otherwise.  Store it as a column in the all_factors dataframe.\n",
    "* Add another similar column to indicate when the month is December"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors['is_January'] = (all_factors.index.get_level_values(0).month == 1).astype(int)\n",
    "all_factors['is_December'] = (all_factors.index.get_level_values(0).month == 12).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weekday, quarter\n",
    "* add columns to the all_factors dataframe that specify the weekday, quarter and year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors['weekday'] = all_factors.index.get_level_values(0).weekday\n",
    "all_factors['quarter'] = all_factors.index.get_level_values(0).quarter\n",
    "all_factors['year'] = all_factors.index.get_level_values(0).year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start and end-of features\n",
    "\n",
    "* The start and end of the week, month, and quarter may have structural differences in trading activity.\n",
    "* [Pandas.date_range](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.date_range.html) takes the start_date, end_date, and frequency.\n",
    "* The [frequency](http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases) for end of month is `BM`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors['month_end'] = all_factors.index.get_level_values(0).isin(pd.date_range(start=factor_start_date, end=universe_end_date, freq='BM')).astype(int)\n",
    "all_factors['month_start'] = all_factors.index.get_level_values(0).isin(pd.date_range(start=factor_start_date, end=universe_end_date, freq='BMS')).astype(int)\n",
    "all_factors['qtr_end'] = all_factors.index.get_level_values(0).isin(pd.date_range(start=factor_start_date, end=universe_end_date, freq='BQ')).astype(int)\n",
    "all_factors['qtr_start'] = all_factors.index.get_level_values(0).isin(pd.date_range(start=factor_start_date, end=universe_end_date, freq='BQS')).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = list(all_factors.columns)\n",
    "features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make Target\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipeline_target = Pipeline(screen=universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "We'll convert returns into 5-quantiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "return_5d_5q = Returns(window_length=5, mask=universe).quantiles(5)\n",
    "return_5d_5q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipeline_target.add(return_5d_5q, 'return_5d_5q')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_df = engine.run_pipeline(pipeline_target, factor_start_date, universe_end_date)\n",
    "targets_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_label = 'return_5d_5q'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors.index.get_level_values(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_df.index.get_level_values(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split into training, validation and test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_into_sets(data, set_sizes):\n",
    "    assert np.sum(set_sizes) == 1\n",
    "    \n",
    "    last_i = 0\n",
    "    sets = []\n",
    "    for set_size in set_sizes:\n",
    "        set_n = int(len(data) * set_size)\n",
    "        sets.append(data[last_i:last_i + set_n])\n",
    "        last_i = last_i + set_n\n",
    "        \n",
    "    return sets\n",
    "\n",
    "def split_by_index(df, index_level, sets):\n",
    "    set_indicies = split_into_sets(df.index.levels[index_level], sets)\n",
    "    \n",
    "    return [df.loc[indicies[0]:indicies[-1]] for indicies in set_indicies]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# put the features and target into one dataframe before \n",
    "# running dropna, so that the rows match.\n",
    "tmp = all_factors.copy()\n",
    "tmp [target_label] = targets_df[target_label]\n",
    "tmp = tmp.dropna()\n",
    "X = tmp[features]\n",
    "y = tmp[target_label]\n",
    "\n",
    "X_train, X_valid, X_test = split_by_index(X, 0, [0.6, 0.2, 0.2])\n",
    "y_train, y_valid, y_test = split_by_index(y, 0, [0.6, 0.2, 0.2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit a random forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import ExtraTreesClassifier, RandomForestClassifier, ExtraTreesRegressor, RandomForestRegressor\n",
    "from sklearn.metrics import log_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = RandomForestClassifier(\n",
    "        n_estimators=10,\n",
    "        max_features='sqrt',\n",
    "        min_samples_split=5000,\n",
    "        bootstrap=True,\n",
    "        oob_score=True,\n",
    "        n_jobs=-1,\n",
    "        criterion='entropy',\n",
    "        verbose=0,\n",
    "        random_state=0\n",
    "    )\n",
    "clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rank features by Feature importance (sklearn)\n",
    "\n",
    "We'll define a function that uses the built in sklearn feature importances, and sorts the features by their feature importance.\n",
    "\n",
    "Note that [numpy.argsort](https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html) returns a list of the original index locations of a list, in the order that would make them sorted in ascending order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argsort([10,30,20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to reverse the order of a list or array is to use the notation `[::-1]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = [3,2,1]\n",
    "tmp[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_importances(m, features):\n",
    "    # TODO: get the feature importances from the model\n",
    "    importances = m.feature_importances_\n",
    "    \n",
    "    # TODO: sort the importances in descending order, and store the indices of that sort\n",
    "    indices = np.argsort(importances)[::-1]\n",
    "    \"\"\"\n",
    "    Iterate through the features, starting with the ones with the highest feature importances\n",
    "    \"\"\"\n",
    "    features_ranked = []\n",
    "    for f in range(X_train.shape[1]):\n",
    "        print(\"%d. %s (%d) (%f)\" % (f+1,features[indices[f]], indices[f], importances[indices[f]]))\n",
    "        features_ranked.append(features[indices[f]])\n",
    "\n",
    "    return features_ranked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### See ranking of features according to sklearn's feature_importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_skl = model_importances(clf, features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Shap library\n",
    "\n",
    "We'll also use the Shap library to determine feature importance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap\n",
    "shap.initjs() #initialize javascript to enable visualizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Shap outputs \n",
    "https://shap.readthedocs.io/en/latest/\n",
    "\n",
    "```\n",
    "shap_values(X, y=None, tree_limit=-1, approximate=False)\n",
    "\n",
    "X:\n",
    "A matrix of samples (# samples x # features) on which to explain the model’s output.\n",
    "\n",
    "tree_limit:\n",
    "Limit the number of trees used by the model.\n",
    "\n",
    "approximate:\n",
    "Run fast, but only roughly approximate the Tree SHAP values\n",
    "```\n",
    "\n",
    ">**For models with a single output this returns a matrix of SHAP values (# samples x # features).** \n",
    "\n",
    ">Each row sums to the difference between the model output for that sample and the expected value of the model output (which is stored in the expected_value attribute of the explainer when it is constant). \n",
    "\n",
    ">**For models with vector outputs this returns a list of such matrices, one for each output.**\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this will take a few seconds to run\n",
    "explainer = shap.TreeExplainer(clf)\n",
    "shap_values = explainer.shap_values(X_train, tree_limit=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The classifier has 5 distinct outputs (5 quantiles).  \n",
    "So the shap_values is a list of 5 numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(shap_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each element of shap_values has one row for each training data point, and one column for each feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# features, data points\n",
    "shap_values[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## local to global feature importance\n",
    "\n",
    "SHAP calculates local feature importance for every training observation (every row).  \n",
    "To calculate global feature importance, take the absolute values of the local feature importances and then take the average across all samples.\n",
    "\n",
    "$GlobalImportance_{i} = \\frac{1}{N}\\sum_{j=1}^{N} |LocalImportance_{i,j}|$  \n",
    "\n",
    "Where there are N samples, and $i$ denotes a particular feature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the built in function to plot the features sorted by global feature importance.  This is taking the average of the absolute values of the shapley values for each feature, to get the global feature importance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean of absolute values for each feature to get global feature importance\n",
    "shap.summary_plot(shap_values, X_train, plot_type=\"bar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the plot shows the first 20 features.  We can write our own function to calculate global feature importance, so that we can see the global feature importance of all features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rank features using SHAP\n",
    "\n",
    "There are a couple classes, one for each quantile.  So the list returned by `shap.shap.TreeExplainer.shap_values()` has one element for each of those classes.  We'll explore how to get the absolute values and then average of those absolute values, for each of the features.  Then we can put this into a function.\n",
    "\n",
    "Remember, here's the formula to aggregate local feature importances into global feature importance:\n",
    "\n",
    "$GlobalImportance_{i} = \\frac{1}{N}\\sum_{j=1}^{N} |LocalImportance_{i,j}|$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can concatenate the 2D arrays in the list `shap_values`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp1 = np.concatenate(shap_values)\n",
    "tmp1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take the absolute values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp2 = np.abs(tmp1)\n",
    "tmp2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take the average for each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp3 = np.nanmean(tmp2,axis=0)\n",
    "tmp3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Implement the function that calculates global feature importance using shapley values, and sorts the features by importance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Challenge: try implementing the function yourself! \"\"\"\n",
    "\n",
    "def model_shap_importances(model, features,X):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also use the starter code below, if you prefer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_shap_importances(model, features,X):\n",
    "    \"\"\"\n",
    "    Note that the observations should be numeric (integer or float).\n",
    "    So booleans should be converted to 1 (True) and 0 (False) \n",
    "    \"\"\"\n",
    "    # TODO: calculate shap values\n",
    "    shap_values = shap.TreeExplainer(model).shap_values(X, tree_limit=5)\n",
    "    \n",
    "    # TODO: concatenate the shap values into one matrix\n",
    "    shap_values_matrix = np.concatenate(shap_values)\n",
    "    \n",
    "    # TODO: take the absolute values\n",
    "    shap_abs = np.abs(shap_values_matrix)\n",
    "    \n",
    "    # TODO: Take the average for each feature (each column)\n",
    "    global_importances = np.nanmean(shap_abs, axis=0)\n",
    "        \n",
    "    # TODO: get the indices sorted in descending order of global feature importance\n",
    "    indices = np.argsort(global_importances)[::-1]\n",
    "    features_ranked = []\n",
    "    for f in range(X.shape[1]):\n",
    "        print(\"%d. %s (%d) (%f)\" % (f+1,features[indices[f]], indices[f], global_importances[indices[f]]))\n",
    "        features_ranked.append(features[indices[f]])\n",
    "        \n",
    "    return features_ranked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this will take a few seconds to run\n",
    "features_ranked = model_shap_importances(clf,features,X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_ranked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion on sector\n",
    "\n",
    "- Random forests can still work with categorical features that are numbers.  For instance, to filter features by sector '5', it's possible for a tree to split on sector < 6, and then sector > 4.  One of the reasons tree-based models are great is because they can still try to interpret data that hasn't been fully cleaned or processed.  However it's still a best practice to one-hot encode categorical features, as this will help to reduce noise, and hopefully help the model's performance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sector category names\n",
    "\n",
    "You'll one-hot encode the sector with sector labels in the project.  Please see some code that can use when you assign category labels to each sector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sector Labels\n",
    "sector_names = pd.read_csv('sector_names.csv')\n",
    "sector_names = sector_names[['Sector','sector_id']]\n",
    "sector_names = sector_names.drop_duplicates()\n",
    "sector_names = sector_names.append(pd.DataFrame([['no sector assigned',-1]], columns = sector_names.columns))\n",
    "sector_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use this dataframe to get the sector name by sector id\n",
    "# here's an example\n",
    "tmp = sector_names.loc[sector_names['sector_id'] == 9]['Sector'].values[0]\n",
    "tmp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-hot encode other features?\n",
    "Are there other features here that you could also one-hot encode?  You can one-hot encode these other features in the project!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## date features\n",
    "\n",
    "The low frequency date parts (end of month, end of quarter etc.) have low importance. This should not be surprising since we are looking at only a couple years of training history and there are only, say, 4 quarters in a year. \n",
    "\n",
    "End of month trading activity may occur some days before the last business day of the month.  To better capture what we think of as end-of the month trading, we can try including the last 5 business days of the month.  Similarly, we can try including the last two weeks of each quarter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Date features helper code\n",
    "You may find some of these functions useful in the project!\n",
    "\n",
    "* We can use [BDay](https://pandas.pydata.org/pandas-docs/stable/timeseries.html) to offset our date_range by a specified number of business days\n",
    "* Also, check out a list of [frequencies](http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases) to choose from, such as `M`, `Q`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas.tseries.offsets import BDay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = all_factors.index.get_level_values(0)\n",
    "tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_1 = tmp + BDay(-1)\n",
    "tmp_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how adding `Bday(-1)` to the DateTimeIndex `tmp` made another DateTimeIndex with the second to last business day of each month."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Union\n",
    "DatetimeIndex has a [union](https://pandas.pydata.org/pandas-docs/version/0.21/generated/pandas.DatetimeIndex.union.html) function that merges one DatetimeIndex with another"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_1.union(tmp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may find this code useful in the project!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prune features\n",
    "\n",
    "In the project, you'll try to improve some of your features, and then check how their feature importance changes (if at all).  Then you can prune your feature list and choose the ones that you want to use in your model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "[solution notebook](feature_importance_solution.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
